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ABSTRACT 

We present a model for plasma heating produced by time-dependent, spatially lo- 
calized reconnection within a flare current sheet separating skewed magnetic fields. 
The reconnection creates flux tubes of new connectivity which subsequently retract at 
Alfvenic speeds from the reconnection site. Heating occurs in gas-dynamic shocks which 
develop inside these tubes. Here we present generalized thin flux tube equations for the 
dynamics of reconnected flux tubes, including pressure-driven parallel dynamics as well 
as temperature dependent, anisotropic viscosity and thermal conductivity. The evo- 
lution of tubes embedded in a uniform, skewed magnetic field, following reconnection 
in a patch, is studied through numerical solutions of these equations, for solar coro- 
nal conditions. Even though viscosity and thermal conductivity are negligible in the 
quiet solar corona, the strong gas-dynamic shocks generated by compressing plasma 
inside reconnected flux tubes generate large velocity and temperature gradients along 
the tube, rendering the diffusive processes dominant. They determine the thickness of 
the shock that evolves up to a steady-state value, although this condition may not be 
reached in the short times involved in a flare. For realistic solar coronal parameters, this 
steady-state shock thickness might be as long as the entire flux tube. For strong shocks 
at low Prandtl numbers, typical of the solar corona, the gas-dynamic shock consists 
of an isothermal sub-shock where all the compression and cooling occur, preceded by 
a thermal front where the temperature increases and most of the heating occurs. We 
estimate the length of each of these sub-regions and the speed of their propagation. 

Subject headings: Magnetic Fields — MHD — Shock Waves — Sun: flares 



Introduction 



Mag netic reconnection has lo n g been accep ted as a significant mech anism for ene r gy release 
in flares (iKopp &: PneumanI 119761 : iPriestI 119811 ). coronal mass ejections (jForbed 12000 : iKlimchuk 



2001 



Low 



2001 



Lin et al 



2OO3I). and as a n interaction process between the solar wind and the 



magnetosphere (jPudovkin &: Semenovlll985l ). It operates in current sheets across which there is a 
discontinuity in the direction of the pre-reconnection field lines. Some kind of diffusion mechanism 
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initiates the connection of field lines from one side of the current sheet to field lines on the other side. 
The exact nature of this process is not yet well understood, but it has recently become clear that it 
must produce an elect ric field restricted to a small portion of the current sheet (jBiskamp &: Schwarz 

Magnetic recqnnectiqn was first modeled in two-dimensions and steady state 



2001: Birn et al.ll2001 



with uniform resistivity (jSweetlll958l : |Parkeiill957l ). but the ener gy release in th is case was too slow 
to explain the rapid development observed in flares and CMEs (lBiskamplll986l ). This problem was 
solved by assuming a 2D steady electric field localized to a small diffusion region with four slow- 
mode shocks (SMSs) attached to it. At these shocks, t he magnetic fie ld is deflected and partially 
annihilated, and th e plasma is heated and accelerated (|Petscheklll964l ). The compressible plasma 
case was studied by ISoward &: Priest] (|l982l ) . 



Localized , non- stead y models, b e gan w ith the work of ISemenov et al.l (jl983l ). followed by 
Biernat et al.l (jl987l ) and iNitta et al.l (120021 ) all of which were two-dimensional. For these spa- 
tially localized (in only one dimension) and short lived episodes, four slow mode shocks (SMSs) 
extend from the reconnection site as well, but they close back together in a finite region that in- 
creases in size as time goes by (as opposed to the infinite length SMSs from steady state). These 
teardrop-sh aped SMSs continue growing and retracting as they sweep up mass, even if reconnec- 
tion ceases ([Semenov et al.lll998l ). Therefore, this conversion of magnetic energy into kinetic and 
thermal energy continues and the dissipated energy can exceed the ohmic dissipation from the 
reconnection electric field, indicating that a small scale, short lived event can affect the global 
structure of magnetic fields. 

There is some observational evidence supporting the basic features of fast magnetic reconnec- 
tion models. The heating and chromospheric eva poration in flare s provides indirect evidence of the 
formation and propagation of slow mode shocks (jLin et al.l 120081 ). Many observational features, as 
plasmoids moving along current sheets can be attributed to unst eady, Petschek-t ype reconnection, 
although they also can be the result of turbulent reconnection (iLin et al.l 120081 ). Current sheets 
themselves are notoriously difficult to observe directly, yet they are a neces sary component of al l 
fast magnetic reconnection models. Those observations providing evidence (jSchettino et al.ll20ld ) 
also suggest high temperatures at the current sheet, consistent with shock heating of material. 

There is also evidence for localized and transient nature of magnetic reconnection. The first 
such evider ice was found in so cal led flux transfer events identified in in situ magnetospheric ob- 
servations (jRussell &: Elphidll978l ). X-ray and EUV observation have recently provided evidence 
of solar reconnection episodes localized in space and time. Dark voids (jMcKenzie Hudsonlll999l ) 
descending from reconnecti on regions are direct ly related to the path of bright loops that appear 
lower down in the corona (jSheeley et al.l 12004 ). These downward-moving features are believed 
to represent flux tub es, created by patchy reco nnection, retracting downward to later form the 
post-eruption arcade (jMcKenzie &: Savagd l2009l ) . 



All of the models above considered reconnection between magnetic fleld lines assumed to be 
perfectly antiparallel. A departure from this special case, called reconnection between skewed fields. 
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introduces a magnetic field component in the ignorable direction, called a guide-field; such models 
are sometimes called 2.5-dimens ional. In this case i t is necessary to supplement the SMS s with 
rotational discontinuities (RDs. iPetschek &: Thornd 119671 : ISowardI Il982l : iHevn et al.lll988l ). The 
effect of the RDs is to produce bulk flow in the ignorable direction in order that the velocity 
change, magnetic fields and shock normal may all lie in a single plane; this is a requirement of the 
SMS often called "co-planarity" . This new flow component is confined to the region between the 
RD and SMS. It is the RDs rather than the SMSs which accelerate the reconnection flux to the 
Alfven speed in an outflow jet. 



Non-st eady reconnec t ion b etween skewed magnetic fields was a nalyzed by ISemenov et al 



(jl992l ) and iBiernat et al.l (jl998l ) for dayside magnetospheric events. iLinton k, Longcopd (j2006l ) 
studied the relaxation of a solar-like low plasma-/3 reconnected flux tube via three-dimensional 
magnetohydro dynamic simulations. Even though the initial configuration is 2.5-dimensional, re- 
connection is initiated within a small patch, producing a distinct bundle of reconnected field lines, 
a flux tube; the ensuing reconnection dynamics is purely three-dimensional. They found that the 
perpe ndicular dyna mics of the reconnected tube was well described by equations for thin flux tubes 
fTFT. ISDruitlll98lh . The equations differed from those used traditionally in convection zone mod- 
els, by their application to the low-/3 corona. Linton and Longcope's study shows RDs propagating 
along the legs of the tube that differ from the usual Petschek SMSs in the sense that they rotate 
the magnetic field without changing its magnitude. In this respect, t hey were like the RDs in 



steady-state 2.5-dimensional models (ISoward 



1982 



Skender et al.l 120031 ). Rather than a coherent 



outflow jet, the Alfven-speed retraction of reconnected flux occurs as a single, transient event. 

Common to all models of fast magnetic reconnection, except those restricted to perfectly anti- 
parallel fi elds, is a second sho ck where plasrn a heat ing occurs, located inside the RDs. In the steady 
models of ISowardI (j 19821 ) and ISkender et al.l (j2003l ). this inner shock is the SMS. At low plasma /3 
and skew angles significantly far from anti-parallel, the SMS differs significantly from a switch-off 
shock: it changes the magnetic field strength ver y little and derives m ost of its thermal energy 



from bulk flow parallel to the field lines. Rece ntlv. iLongcope et al. (120091 1 presented a TFT model 
which included flow parallel to the fleld lines ( Linton &: Longcope 20od . had not investigated the 



parallel plasma dynamics). They found a shock of apparently different character, a gas dynamic 
shock (GDS), changing only the parallel flow and leaving the magnetic field strength completely 
unchanged. Upon further consideration, however, this shock appears related to the nearly-parallel 
SMS in 2.5 steady models. In fact the post-GDS temperatures in the TFT model are almost 
identical to the post-SMS temperatures in steady-state MHD models. In both kinds of models, 
plasma heating occurs through a two-stage process notably different from commonly-used switch- 
off models, applicable to purely anti-parallel reconnection. First the RD converts magnetic energy 
to bulk kinetic energy in the form of Alfven-speed fiow partly parallel to the field line. The parallel 
component of this flow, present in both, is then thermalized in a second shock of large hydrodynamic 
Mach number. 



Also common to all the reconnection models is the presence of Alfvenic outflows originating 
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at the reconnection site. In steady state models these f orm coherent j ets, for which there is some 
observational support, although it is typically indirect (jTsunetal 119971 ). In transient reconnection 
models any outflow is more accurately described as a distinct retraction event rather than a steady 
jet. Observational evidence fo r these events is provided by the ongoing investigations of downfiowing 
featu res above flare arcades ([McKenzie &: HudsonI Il999l : ISheelev et al.l 12004 : iMcKenzie Sz Savage 
20091 ). While few if any of these events have measured speeds matching the local Alfven speed, 
they are generally interpreted as evidence of reconnection. Either these outflows/downflows were 
Alfvenic at some earlier time or all reconnection theories must be revised. In the model presented 
here we continue to assume reconnection outflow at the Alfven speed. 

While all of the foregoing reconnection models predict shocks, few investigations of them have 
explicitly included viscosity or thermal conductivity. This omission is noteworthy since strictly 
speaking the irreversible nature of any shock demands a non-ideal effect; in the collisional plasma 
of the solar corona this is usually assumed to be viscosity. Nevertheless reconnection models gener- 
ally neglect viscosity altoge ther to say nothing of a ccounting for its s trong temperature dependence 
and directional anisotropy. ICraig et al.l (j2005l ) and iLitvinenkd (j2005l ) studied the effects of tem per- 
ature independent v i scosit y in various reconnection disturbances. Only recently, ICraiei (|2008l ) and 
Craig Litvinenkd (j2009l ) included anisotropic and temperature dependent viscous dissipation in 
a three-dimensional steady-state model at X-points. 

For strongly magnetized plasma, viscous diffusion and thermal conduction are highly anisotropic. 
They are directed along the magnetic field since perpendicular coefficient are many orders of magni- 
tude smaller than their parallel counterparts. Transport coefficients in a fully-ionized plasma each 
have a strong dependence on temperature (~ Ta) making transport far more significant at high 
post-shock temperatures than it is ass umed to be in general. It has therefore been recognized as a 



potentially significant factor in flares (| Car gill et al 



1995 



Forbes &: Malherbelll986l ). In particular 



the heat generated at SMSs is predicted to be conducted along rec onnected field lines ahead of other 
reconnection effects, possibly driving chromospheric evaporation ([Forbes et al.lll989l ). Yokoyama &; 
Shi bata (1997, 2001) included anisotropic thermal conduction in 2D simulations of reconnection, 
and lSeaton &: ForbesI (|2009l ) presented a detailed theoretical analysis of these simulations. In most 
of these two-dimensional, anti-parallel models, the heat conduction front extends in front of SMSs 
as expected. 

Here we seek to understand the role played by anisotropic transport in transient (non-steady) 
reconnection between skew ed n iagnetic fields. Toward this end we generalize the TFT equations of 
Linton Longcopd (|2006l ) and lLongcope et al.l (j2009l ) to include pressure-driven parallel dynamics 
as well as temperature dependent, anisotropic viscosity and t hermal conductiv i ty. T he latter are 
essential for self-consistently producing the shocks predicted by iLongcope et al.l (120091 ) analytically 
using conservation laws. 

After presenting the general equations, we consider the flare-like scenario depicted in Fig. [1] 
where reconnection occurs within a current sheet above an arcade of post-flare loops. Short lived. 
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patchy (fast) reconnection (depicted in the Figure as a small sphere) is assumed to reconnect the 
skewed field lines from opposite sides of the sheet. The small patch creates a small bundle of field 
lines that subsequently retract at Alfvenic speeds. We approximate the magnitude of the magnetic 
field as uniform and assume strongly magnetized plasma (low plasma-/?). With these assumptions, 
modified thin flux tube equations are presented and the evolution of a thin flux tube that has been 
reconnected is analyzed in detail. 



We solve our new equations numerically for the same reconnection scenario studied bv lLongcope et al 



([20091), using realistic solar coronal parameters. The diffusive processes determine the thickness and 
in ternal structure of the shock s. We analyze this in t ernal struc ture using met hods similar to those 
of lThomasI |l94J) , ICradl ^195^ ) , Icilbarg k Paoluccil , and Ixennell (|l988l l . Our own analysis is 

novel in its focus on solar flare conditions, including the very low Prandtl number (ratio of viscosity 
to thermal conductivity) and strong temperature dependence. In particular we find that for strong 
shocks at low Prandtl number, the gas dynamic shock comprises a narrow iso-thermal sub-shock 
behind a very large conduction front. Due to the strong temperature dependence of the coefficients 
the conduction front has significant extent; far greater than the term "internal structure" might 
suggest. It is related to the conduction front anticipated in two-dimensional anti-parallel models, 
but here we offer a detailed analysis of its structure. 

In addition to thermal conductivity we find it it is necessary to include viscosity in order to 
account for the isothermal sub-shock region. Most of the pl asma co r npres sion occurs isothermally 
within the sub-shock so the plasma is actually cooled by it. iFowlesI (| 19751 ) showed that for super- 
cri tical shocks, it i s impossible to assume inviscid behavior, no matter how large thermal conduction 



IS. 



Coronitil (jl970l ) studied dissipation discontinuities in shock waves for temperature independent 



transport coefficients. The seemingly harmless assumption of temperature independent coefficients 
leads to shock thicknesses on the order of the particle's mean-free-path, raising concerns about the 
adequacy of a fluid treatment: when th e thickness of th e shock is of the order of particle mean 
free-path, the fluid model breaks down (jCampbelll 11984 ). When transport coefficients are given 
an accurate temperature dependence, however, we find the thickness of the shocks exceeds the 
particle mean free paths by at least an order of magnitude. Our simulations also show that the 
ratio between the heat flux and the nominal free-streaming heat flux is smaller than 0.004. These 
two results allow us to treat our plasma as a collisional fluid. 

The layout of the present paper is the following. The next section is a brief review of the 
general transport effects in magnetohydrodynamics. The following section describes the thin flux 
tube assumptions that lead to the final thin flux tube equations presented in the third section. 
Section H] presents simulations of the retraction of a thin flux tube reconnected in a localized region 
inside a current sheet with uniform, skewed magnetic flelds for a case with low plasma-/3, and 
parameters relevant to the solar corona. In the fifth section, an analytical analysis of the inner 
structure of the shocks is presented. Finally, in the last section, the main results of the paper and 
possible observational implications of the inner structure of the shocks are discussed. 
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2. MHD equations 



When resistivity, radiation and gravity are neglected, the magneto-hydrodynamic (MHD) equa- 
tions fora charge-neutral plasma (we will assume sufficient collisionality to justify their use) are 



(|Shulll99ll . chap. 22): 



1 Dp 
Dv 

'dt 



-V-v, 

-V(P + g)+V. 



BB 

47r 



v-n, 



Dei PDp ^ ^ 

p — ^ = ^'-V-q, 

^ Dt p Dt ^' 

-V X (B X v) . 



(2-1) 
(2-2) 
(2-3) 
(2-4) 



Here, p, P, B, v, T, s, ej, H, and q are the density, plasma pressure, magnetic field, velocity, 
temperature, specific entropy, specific internal energy, viscous stress tensor, viscous heat, and 
thermal conduction flux of the fluid, respectivel y. The adyectiv e derivative is defined as D /Dt = 
d/dt -|- V • V. In Eq. (I2-2|) . dyadic notation (iGoldsteinI Il959l . chap. 5) is used for the second 
right-hand-term. 

When the plasma is strong ly magnetized by a field taken in the z direction, the stress tensor 
takes a simple gyrotropic form (jBraginskiil 1 1 965l ) 



where the zz component of the rate-of-stress tensor is 

^dvz 



dz 



3 



(2-5) 



(2-6) 



Here, r/ is the largest of the dynamic viscosity coefficients, determined essentially by the ions (the 
electron-to-ion viscosity ratio is ?ye/f?i — 0.01). The strongly magnetized limit of the conductive 
heat fiux is 



(2-7) 



where Vy = zz • V is the derivative parallel to the magnetic field, and k is the parallel thermal 
conduction coefficient, determined mostly by the electrons (the electron-to-ion heat conductivity 
ratio, Ke/ni ~ 24). 



We close the system using the ideal gas law 

P : 



^ = 2nfc^r, 
m 

P 



(7-l)p' 
P' 



cv In 



+ So, 



(2-8) 
(2-9) 

(2-10) 
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where 7 is the adiabatic gas constant, rn the average particle mass (for a fully ionized Hydrogen 
plasma, it is approximately equal to half the mass of the proton), n is the electron number density, 
cv = / {"f — is the specific heat at constant volume, and sq is a constant. 



3. Thin Flux Tube Equations 

3.1. The Thin Flux Tube Assumptions 

We consider an untwisted, thin tube of magnetic flux </>, like the one depicted in Fig. O 
embedded in an strongly magnetized external fluid. The external magnetic field, thermal pressure 
and density will be denoted Be, Pe, Pe, respectively. These quantities may depend on position, 
but not time, since we will assume the background configuration to be in equilibrium. The tube is 
considered so thin that the background quantities remain unperturbed as the tube moves through 
the external plasma. It is also tine compared to the local radius of curvature. 

We will assume that the tube can be represented by an internal field line, R(Z), the axis (see 
Fig. [2]). The field is untwisted in the sense that all field lines are parallel to the axis. The time 
it takes for a fast magneto-sonic wave to establish pressure balance with its surroundings scales 
as = a(/)/y^c^(/) + v'^{l), where Va is the Alfven speed, and Cg is the sound speed. Since the 
fast-magnetosonic time is negligible with respect to the Alfven time over a characteristic scale L^h, 
TA = Lch/va{l), we assume instantaneous, local pressure balance 

P+— = Pe + —, (3-1) 

where, P and B are pressure and magnetic field inside the tube. 

It is convenient to define a Frenet coordinate about the axis R, parameterized by its arc-length 
I, as shown in Fig. [2j The Frenet unit vectors defining the coordinate system are 

^^(^) - mm - 

m = k(/)re(/) = ^r,(/), (3-3) 
a(/) = p(/) X b(/). (3-4) 

The unit vector b is parallel to the axis, p is the normal vector and n is the bi-normal completing 
the system. Local derivative operators associated with the unit vectors are 

5|| = b-V = 9;, 
V_L = -b X (b X V). 



We assume the flux tube is sufficiently isolated from its surrounding that no viscous stress 
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crosses its surface. The stress tensor is therefore 

n = -v(h-^^hh. (3-5) 

For parameters typical at the high corona, T = 1 MK, n = 10^ cm"^, B = 10 G, L = 100 Mm, the 
Reynold's number is i?^ ~ 3 x 10'^, indicating a relative small contribution from viscosity. For such 
high Reynolds numbers, the viscous stress tensor may become sig nificant only at steep velocity 
gradients, like shocks. Shocks in thin flux tubes were predicted by ()Longcope et al.ll2009l ). 



The small Prandtl number of an ionized plasma, Pr = rj/nm ~ 0.01, indicates that the effect 
of thermal conduction is larger than viscosity. Temperature gradients will then be even smaller at 
shocks than velocity gradients. 



3.2. Thin Flux Tube Equations 

The thin flux tube equations are derived by applying the MHD equations to a small piece of 
the tube with mass 5m., like the one shown in Fig. [2j The position of every point in the tube can 
be expressed as r' = R + r. Vector r belongs to the circular area perpendicular to the tube axis, 
and can be expressed as r = r sin n + r cos p = rr, where r and 9 are the polar coordinates 
associated to the n-p plane. 

The local radius of the tube, a(/), is small compared with scales of variation along the tube, 
permiting Taylor-expand any state variable / as follows 



/(r) 



E 

fc=0 



[(r • V) V] (0 



k\ 



fc=0 j=o 



cos I 



Xsin( 



(3-6) 



The second equality was obtained expanding the binomial term in the local polar coordinates. 
Thinness also allows us to assume that the piece's axis lays on a plane, and its radius of curvature 
is constant along the small piece. For this case, the differentials of area and volume for the piece 
can be expressed as 



dA 



dAi 



dv 



h{l) r dr dO, 



dr' dr' 
'de^'dl 



?(/) 



ail) 



dr' dr' 
'de^'dl 
sin( 



dl de, 

r=a{l) 

a^il)sme 
dr' 

• — dldO dr 
dr 



-h{l)a{l)^^^dld9, 



dl dO dr, 



(3-7) 



(3-8) 



(3-9) 
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An example of the expansion (j3-6p is the mass of the tube piece. Computing the volume 
integral using expression (j3-9p for the volume differential, and expanding the density around the 
center of the tube, gives 

6m = [ p dv 

/u\ 1 1-1+51 ak„ri/\ rail) r2-K 

^k+l ^ 



k=0 i=0 



JO 



X {cos0)''-'{sm9y f 1 _LJ^] dOdr. 

1+51 A+Sl 

p{l')TTa\l') dl' + / 0{a\l'))dl\ 
I Ji 

where 0{a'^{l')) represents an error of fourth order in the local radius. If the arc-length interval, 
61, is small then the above equation can be expressed as 

6m{l) ~ p{l)TTa^{l)6l, (3-10) 

Expanding the dot product B • b, gives a similar expression for the magnetic flux 

B • dAii 

B{l)Tra\l) + 0{a\l)) ~ B{l)Tra\l). (3-11) 



For a given tube piece, the ratio between the mass and the magnetic flux, 6p = 6m/ (p = 
{p/B)6l, is constant. This definition, along with Eq. (|3-10p and (|3-lip . allow us to define the 
integrated mass per flux, p = J {p/B)dl, to parameterize the tube instead of I. The differential 
with respect to p is 

d B d , , 

a-, = 7a- 

In order to obtain the thin flux tube mass equation we integrate over the piece's volume every 
term of the corresponding MHD equation, and keep only the first non-vanishing terms. To first 
order, the velocity at any point in the tube is the sum of the velocity of the center of the tube plus 
a radial term Vr- This radial velocity at the surface of the tube is equal to the change in local radius 
[D / Dt)a{l). Finally, applying the fundamental theorem of calculus to recover integrals along the 
arc- length, and expression (j3-lip . the final result is 

D (B\ ^ dv , , 



In the integration of the momentum equation, (j2-2p . the surface integrals for the perpendicular 
area vanish since B • dA_|_ = and 11 • dA_|_ = (in this last equality we have used definition (j3-5p ). 



- 10 - 



For the first term of the right hand side, constraint ()3-ip is apphed and expanded around the center 
of the tube yielding 



dP 



+ 2k Pe - P + 



Svr 



+ 



B d 



b?7 5v 
"b V '~dl 



(3-14) 



where every term is evaluated at the center of the tube. In the above equation, the viscous 
momentum term is a perfect derivative in the integrated mass. There is no viscous momentum 
interchange between the tube and its surroundings. 

To determine the form of the viscous heating, we compute the rate of change in kinetic energy 
as the dot product between the momentum Equation (|3-14|) and the velocity. Balancing this loss 
with a heating rate gives 



V b 



'dl 



(3-15) 



With the above definition, and similar treatment than for the mass and momentum equation, the 
thin flux tube energy equation becomes 



Dei 
15t 



P Dp r] 



b 



5v 

In 



+ 



B d / K d 



p dl \Bdl 



(3-16) 



where we see that the thermal conduction term is also a perfect derivative in the integrated mass 
mass per flux (there is no thermal conduction across the tube). 

For our simulations, it is more convenient to re- write the energy equation by deflning a new 
variable Pc in the following way 



P{p,t) 



p{p,0) 



Peip) 



(3-17) 



In an adiabatic case, Pc is constant in time. The change in plasma entropy with respect to the initial 
state (assumed to be in equilibrium) is directly related to Pc, where As = cvl'n.[Pcip,t)/Pc{p)]. 
With definition (I3-I7p . the energy equation is transformed to an equation for Pc 



DPc 
Dt 



(7 



1)1^ 



(3-18) 



When density is known, the above equation allows us to update pressure, and use it in the momen- 
tum equation. 
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Heating and cooling occur where there is a change in entropy. From Eq. (|2-3p . ()3-16p . and 
(j3-18p . the volumetric heating rate can be written as 



Ds 1 pY DPc 



Dt 7 - 1 VPe/ Dt 



In this equation, the viscosity term is always positive, and therefore always heats the plasma. On 
the other hand, the thermal conduction term can have either sign, depending on the temperature's 
arc-length second derivative. 



The transport coefficients in an ionized plasma depend om temperature as (|Spitzerlll962l ) 



V = VcT^, (3-20) 



5 



K = KcT2, (3-21) 

where rjc and Hc are constants. 

For thin flux tubes, there is no need to solve the induction equation, since the magnitude of 
the magnetic field inside the tube is determined by constraint ()3-ip . and its direction is given by 
the position of the tube piece governed by the momentum equation (resistivity has been neglected, 
and therefore the magnetic field follows the fluid). 

It is possible to describe the TFT equations in dimensionless variables, by defining a charac- 
teristic magnetic field magnitude, Bq, density, po, and length, Lq. For example, these characteristic 
variables may represent the magnetic field magnitude and density at the reconnection point, x/j, 
and the length of a typical reconnected tube, respectively. Then, Bq = Bi.{xji), and po = pe(xif,). 
With these values we can construct an associated Alfven speed v^q = Bq/ y/Anpo, and an associated 
temperature = v^Qfn/kB, to non-dimensionalized velocities and temperatures. 

The new dimensionless variables are 

,,* -Bp * _£_ J.* 4 VAO J* p* p 47r p* p 47r rp* T 

v* = ^, B*=B^, x* = ^, k*=kLo, b*=b, 
^* = OT^ = (T*)^/^with = ,.gg^. 



J" = K* (r*)^/2 with K* = ^ /"(y)^^' . 



(3-22) 



With these new variables, the dimensionless thin flux tube equations have exactly the same form 
as Eq. (j3-13p . (|3-14p . ()3-16p . and (j3-18p . except that the Boltzmann constant is no longer present. 
From now on, we will refer as the thin flux tube equations as the dimensionless ones. 
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4. Solutions of Reconnection Dynamics 

4.1. Low Plasma-/? Uniform Background Case 

The results of the last section are quite general. In this section, we will simplify the equations 
to apply them to the low-/? solar corona. For simplicity, we will also assume uniform background 
magnetic field. 

We will consider skewed, uniform external magnetic field forming an angle 2^ on opposite sides 
of a static, infinite current sheet, as shown in Fig. [3l The current sheet is assumed to be infinitely 
thin, represented by the x-y plane. The initial pressure is also assumed to be uniform everywhere. 
A spatially and temporally localized reconnection event has occurred (by an unspecified physical 
mechanism), connecting some field lines from one side of the current sheet, to some on the other 
side. The amount of fiux is determined by the size of the reconnection region. This bundle of field 
lines forms two V-shaped fiux tubes (as shown in panel (b) of Fig. [3]) and magnetic tension at 
their cusps causes them to retract. These tubes will remain as coherent entities, preserving their 
magnetic fiux 0, while they move. No further reconnection is assumed. 

For this uniform case, the characteristic magnetic field magnitude and density at the reconnec- 
tion point are the same as everywhere else. Then, Bq = Be, po = Pe, and Lq is some characteristic 
length. The perpendicular gradient of the external total pressure is zero. 

In the solar corona, the ratio between thermal pressure and magnetic pressure, /3 = SttP/B'^, is 
small. For a plasma with magnetic field B^ = 10 G, number density rig = 10® cm"^, and temperature 
Te = 1 MK, (3 is approximately equal to 0.007. Consequently, we will neglect the thermal pressure 
with respect to the magnetic pressure. The initially small internal pressure might increase due 
to compressions along the tube, but we will consider it small compared to the internal magnetic 
pressure at all times. In the parallel direction, there is no force related to the magnetic field, 
therefore the parallel derivative of the thermal pressure cannot be neglected. 

With the above considerations, and making explicit the temperature dependence of the trans- 
port coefficients, we can re- write the constraints and thin fiux tube equations in the following 
way 
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where 



= + £K + £M, 

with specific internal, kinetic, and magnetic energies 

P 



(7-l)p' 



|v| 



47rp 



(4-5) 



(4-6) 



The specific magnetic energy of the tube doubles because there is work done by the background 
magnetic field expanding into volume vacated by motion of the tube. We have assumed that the 
external field returns to its original value after the tube has passed. In more realistic scenarios this 
value may be slightly different. 

Conservation of energy improves stability of any numerical solution. Toward this end, we will 
neglect the small last term in Eq. (l4-4p (proportional to /?), that is not part of the perfect derivative 
in the integrated mass. With this simplification, Eq. ()3-18|) remains unchanged, and the change of 
the tube's total energy per flux is equal to 
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(4-7) 



where /Uq ^ind represent the end points of the tube. If the end points are fixed, and there is no 
temperature gradient there, the total energy of the tube is conserved. The neglected term came 
from the change in kinetic energy, calculated from Eq. (|4-3p . and to be consistent, the correspondent 
perpendicular term will be neglected in the momentum equation, as well. Then, the change in total 
momentum per flux of the tube is equal to 
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(4-8) 



With the neglect of the proportional to /3 small term, and the fa ct that Pp is uniform, t he ideal part 
of Eq. (|4-3p reduces to the thin flux tube equation presented by lLongcope et al.l (j2009l ). Therefore, 
gas-dynamics shocks and rotational discontinuities are expected along the reconnected flux tubes. 



The volumetric heating rate when the magnetic fleld magnitude is uniform becomes 



Q 



"11 



2 s^ri 



(4-9) 



In the above equation, viscosity always contributes to heat the plasma. Depending on the curvature 
of temperature to the 7/2 power, and the value of the Prandtl number, the thermal conduction 
term may contribute with plasma cooling or heating. 
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4.2. Simulations 

We developed a computer program, called DEFT (Dynamical Evolution of Flux Tubes), that 
solves the thin flux tube equations for the retraction of the two reconnected thin flux tubes described 
in last section. The program takes into account the effect of the transport coefficients, including 
their strong dependence on temperature. It also can perform simulations using real coronal dimen- 
sionless parameters, as well as implement a strong anisotropy in the transport coefficients (heat 
conduction and momentum transport are allowed only in the direction parallel to the magnetic 
field). This is far more difficult to achieve in fully 3D MHD programs. The DEFT program can 
be implemented for different background conditions as well, but we will restrict ourselves in this 
paper to the low /3 uniform background case from last section. In this case, the two tubes moving 
in opposite directions are symmetric. In cases with non-uniform background that might not be so. 
We use the conservative form of the thin flux tube equations introduced in the last section. 

The DEFT program uses a staggered mesh where each tube piece is represented by grid points 
at its ends, and by its mass per flux. The code implements a Lagrangian approach where each tube 
piece is followed, guaranteeing mass conservation. Due to mirror symmetry, only half of the tube 
is simulated, the other half is a mirror image in the x-direction. 

The initial state corresponds to two opposite V-shaped flux tubes, connected at the reconnec- 
tion region, as shown in Fig. [3l The initial angle between the magnetic field lines at each side of 
the current sheet is equal to 2^. Each tube is divided in N segments separated by A^-|- 1 tube-point 
positions. The initial tubes are clearly out of equilibrium due to the sharp angle between the field 
lines. To avoid introducing length scales at the limit of resolution, the central parts of the initial 
tubes are smoothed. The end-points of the tube are assumed to be fixed, and to have no temper- 
ature gradient with the external boundaries (this last condition ensures no heat transfer from the 
end-points). 

Since the reconnected tubes originated as background plasma they are initialized with uniform 
density, p = pe, and thermal pressure P = Pe- The mass per fiux of each piece is equal to 5m = 
PeSle/Be, where 5le is the segment's initial length. The magnetic field magnitude, B^, is constant 
everywhere throughout the entire simulation. The initial velocities are zero (equilibrium state), even 
for the central bend (we are not interested in the details of the initial transient evolution). Viscosity 
and thermal conductivity are calculated using Eq. ()3-20p and ()3-2ip . where T = P/p (unitless 
quantities). The unit parallel vectors are calculated at segment's centers from definition p-2p . and 
curvature vectors are calculated at grid points from definition (|3-3p . Arc- length derivatives are 
computed using centered differences. 

At every time, the -|- 1 positions and velocities, as well as the A^ pressure pre-factors Pc, 
are advanced using the current state variables in Eq. (14-30 and p-lSp . following a simple Eularian 
approach, for example v(t -|- 5t) = v(t) + (D/Dt)v(t)6t. The time interval is chosen to satisfy the 



Courant-Friedrichs-Lewy condition (jCourant et al.lll967l ). 
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All terms in Eq. (j4-3p are evaluated at grid points, while terms in Eq. ()3-18p are evaluated 
within segments. The state variables and the unit tangent vector b are defined inside the segments. 
They are averaged between neighbors when needed at grid points. A similar approach is used 
when vectors defined at tube points are needed at segment centers. From the advanced positions, 
velocities, and pressure pre- factors Pc, it is possible to define all the advanced state variables. The 
density of each segment at a given time is calculated as p = {5m/ ^) (Be/ 61), where 61 is the updated 
length of the segment. Pressure follows from Eq. (j3-17[) and (13-18p . and the process is repeated 
until the end of the simulation. 

As the two tubes retract they move in opposite directions, but due to the symmetry of the 
background configuration their dynamics are identical. Therefore, we present results only for the 
downward moving one. The simulated equations are dimensionless, but for concreteness the results 
are re-dimensionalized using characteristic magnetic field magnitude, density and length. For this 
illustrative purpose we choose Bq = 10 G, density po = ngm with electron number density no = 10^ 
cm"^, and rn equal to half the mass of the proton. The initial angle between the pre-reconnection 
field lines is 2^ = 90 degrees. The associated Alfven speed and temperature are vao ~ 2.2 Mm/s, 
and ~ 290 MK. We will consider a plasma at Tg = 1 MK, and the plasma- /3 in this case is 
equal to /3 = 2Te/TQ^ ~ 0.007. For a typical length of Lq = 100 Mm, the Reynolds number is 
Rj^ ~ 2982. An Alfven time (time to move one characteristic length at the Alfven speed) would 
be approximately 46 s. As mentioned before, the Prandtl number of the solar corona is typically 
equal to 0.01, and we will use this value in our simulations. The total number of segments in the 
tube is 4000. 

The simulation is carried out for approximately 3200 s, to allow for the full development of 
shocks (discussed in detail in the next section). The initial length of the tube was chosen to be 
unrealistically long in order to achieve the steady state for the shocks. As we will see in the next 
section, the thickness of the shocks, when fully developed, is of the order of the length of an entire 
coronal loop. The fluid inside the tube evolves towards the steady state, which it probably would 
not attain within a realistic flare time. 

The left panel of Fig. S] shows the downward moving tube at different times. The tube is 
seen to be composed of three main segments, one horizontal segment moving downward at the 
y-projection of the Alfven speed, and two legs from the unperturbed parts of the initial tube. The 
two corners where these three pieces come together move at the Alfven speed along the legs of 
the initial tube and will be called bends (triangle shapes in the Figure). In the right panel, the 
y-position of the right bend is plotted as function of time for the same times as in the left panel. 
The solid line correspond to a line with a slope equal to the y-projection of the Alfven speed at the 
reconnection point (vAe)- 

The magnitude of the magnetic field remains constant along the tube due to constraint (|4-ip . 
therefore the area of the tube also remains constant (Eq. |3-llj ). At the bends, the initially 
stationary plasma is deflected along the bisector of the angle between the two adjacent straight 
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segments (average direction of the curvature force), as shown in Fig. [5j In this Figure, a time 
approximately equal to 2550 s was chosen to show the shape and internal structure of the tube long 
after reconnection. 



The velocity of the deflected plasma (region 2 in Fig. [5]) is (see iLongcope et al.ll2009l ) 



V2 = -2vAe sin^ X - VAe sin(C)y. (4-10) 

This velocity is depicted in the top panel of the mentioned Figure (the left side is a mirror image). 

The parallel component (x direction) of this velocity is proportional to the Alfven speed, and 
is therefore generally supersonic {vAe = Cse\/2/7/3). The tube gets shorter as it moves, creating 
parallel supersonic flows at the bends. These supersonic flows are directed toward the center of the 
tube where they collide. The small, but non-zero thermal pressure there prevents plasma from piling 
up, and stops it in the parallel direction. The two colliding flows generate two gas-dynamic shocks 
(the magnetic field remains constant across these shocks) that move outward from the center of the 
tube. In the top panel of Fig. \5\ velocity discontinuities in the parallel direction between region 
1 (post-shock region) and 2 (inflow region), are evident. The vertical ticks indicate the beginning 
and end positions of these shocks. In the bottom panel of the Figure, the left side depicts the 
temperature profile along the tube, and the right side, the plasma-/3 profile, both for t = 2550 s. 
We see that temperature rises almost one order of magnitude from the pre-shock value, and the 
plasma-/3 more than 30 times its initial value, and neither changes at the bends since these corners 
are Alfven waves. 

From Eq. ()4-7p and the boundary conditions at end-points of the tube, we see that the total 
energy of the tube is conserved (as shown in Fig. [6]). Even though, temperature changes significantly 
at the shock, only a small fraction (less than 10 %) of the available magnetic energy is converted to 
thermal energy (Fig. [6|). The magnetic energy of the tube decreases with time as the tube shortens 
(the magnitude of the field remains constant in the tube, but its length decreases), and is mostly 
converted to kinetic energy at the bends. 



5. Gas-Dynamic Shock Internal Structure 



Transport coefficients (viscosity and thermal conduction) and their temperature dependence 
determine the inner structure of the gas-dynamic shocks. As the flux tubes retract, they follow 
a time dependent evolution until they reach a steady-state characterized by a trai isition between 
upstream and downstream states dictated by the Rankine-Hugoniot conditions (jRankine 



187C 



Hugoniotl 118871 ). The diffusive terms smooth out the transitions but do not alter the upstream or 
downstream values provided these regions are asymptotically uniform {dv/dl — )• and dT/dl — )• 0). 
These asymptotic values are governed by conservations laws observed by the diffusion. 



Gas dynamics shocks in reconnected thin flux tubes were introduced bv lLongcope et al 



(120091 ) ■ 
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in the case of ideal thin flux tube equations. It was shown there that only two kind of discontinuities 
are supported by thin, reconnected thin tubes embedded in a low plasma-/3 uniform background: one 
that corresponds to a bend moving at the Alfven speed, where temperature, density and pressure 
are continuous, and a second one that corresponds to a gas-dynamic shock within a straight section 
of the tube. Both of these discontinuities were present in our simulation. In this section, we will 
focus on the later one, where heating and cooling occur (the heating at the bends is negligible). 

To determine the effect of diffusive terms in t he inner structure of the shocks, we follow an 



appr o ach siinilar to the one used by several authors (|Thomaslll944l : lGradlll952l : lGilbarg &: Paolucci 
19531 : iKennellllQSg ). We consider a straight section [dh/dl = 0) and further assume the presence 
of a transition (shock) moving to the right at constant speed Vg- The inflow region will be denoted 
by a "2" subscript, and the post-shock region by a "1" subscript. In our reconnection scenario, the 
inflow region corresponds to the plasma moving supersonically after being accelerated by the bend 
(region 2 in Fig. [5]), and the post-shock region corresponds to the center of the tube (region 1 in 
the same Figure). 



5.1. Steady-state Solution 

The Mach number (the ratio between the fluid velocity and the local sound speed) in the 
i nflow region is deter mined by the initial plasma- /3 and the initial angle between the field lines 



(jLongcope et al.ll2009l ). namely 



M, = ,/ — sin^^l), (5-1) 



and the Mach number in the reference frame of the shock is 

M2,s = M2 + — = J^, (5-2) 

where Cs,2 is the sound speed in region 2, and U2 and T2 are the pre-shock velocity in the shock 
reference frame, and temperature, respectively. For our simulation in the last section, the Mach 
numbers were M2 — 3.85, and M2^s — 5.3. 

In the reference frame of the shock, the inflow velocity corresponds to plasma moving at 
U2 = 'Vx,2 — Vs (< 0), and the post-shock flow is ui = —Vg (plasma is stopped in the parallel 
direction). The ratio between the shock speed, and the inflow speed magnitude, |fx,2|, is given by 



\vx2\ V M4 16 



l^+(l±i)!-i^. (5-3) 



If constant states are assumed, d/dt = for all the quantities calculated in the moving frame of 
this discontinuity. In this reference frame, the steady-state, conservative, unitless TFT equations 
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for uniform background case and low beta can be written as 



piui = P2U2 = pu, (5-4) 
Piuj + Pi = P2ul + P2 = pu^ + P-ri—, (5-5) 



Piui + mPi = p2u-^ + U2P2 



7-1 ' " 7-1 

o „ 27 du dT , , 

= pu^ + uP '--2u7]—-2k—, 5-6 

7 — 1 ol 01 

The second equalities in each of the above equations refer to the state variables inside the shock. 
The usual Rankine-Hugoniot conditions can be found from the first equalities in each of the above 
equations. 

It is convenient for our analysis of the internal structure of the shock, to define a characteristic 
velocity u* = ^ {u2 + T2/U2) = (mf^) (l + l/l^^^l s)- With this definition, the Rankine-Hugoniot 
conditions can be written as 

P2 Ui {AjU* /-f + 1) - U2 , , 

— = — = , (5-7) 

Pi U2 U2 

Ti ui{2u*-ui) 



T, 



2 



(5-8) 



The plasma-/? ratio can be calculated as /3i//32 = {pi/ p2){Ti/T2). In Fig. [5l [H and [TOl horizontal 
dashed lines show the theoretical steady-state jump values for our simulation. Even-though the 
post-shock plasma-/3 is more than 30 times the initial value, it is still below unity. 

The second equalities in Eq. (I5-5P and (l5-6p . are differential equations for the shock's inter- 
nal state variables. Their boundary conditions are given by the steady-state Rankine-Hugoniot 
conditions. Combining these two equations we obtain 

2(7 - 1) dT ^ 2T - u'^jj - 1) + 4(7 - l)u*u + + 1) - 4jU2U* ^ h{T, u) 
Pr u du — 2uu* + T g{T, u) ' 

where we have defined functions h(T, u) and g{T, u) as the numerator and denominator of the above 
differential equation. indicates the Prandtl number. This equation is valid for any temperature 
dependence of the transport coefficients, as long as both are the same; and can be numerically 
integrated for different Prandtl numbers and shock strengths. Numerical solutions of this equation 
in conjunction with Eq. (15-5P or (I5-6P provide a method to determine the thickness of the shock. 

In Fig. [71 numerical solutions to the above equation are plotted for the same Mach number 
as in our simulation from last section (temperature and velocity are normalized to the pre-shock 
values). The pre-shock and post-shock values are fixed by the strength of the shock, independently 
of the transport coefficients, therefore we will call this kind of dia grams fixed po int diagrams (they 



correspond to a v elocity contra st-temperature plane described bv iKennell (jl988l ). or Navier-Stokes 



direction fields by iGradI (119521 )). Each solid line corresponds to a numerical solution to the above 
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equation for a given Prandtl number (increasing thickness of the line represents increasing value in 
the Prandtl number). In this Figure, the mixed dotted and dashed line shows the curve h{T, u) = 
that corresponds to the limit of infinite Prandtl number P^. As the Prandtl number increases 
from low values, the solutions (solid lines) approach this limit case. The dashed line represents the 
g{T, u) = curve, and asterisks indicate grid points (inside the shock region) from our simulation 
{Pr = 0.01) for the same time chosen for Fig. [5l Simulations coincide with the steady-state 
theoretical predictions. There are several grid points in each sub-region of the shock, which shows 
that the DEFT program can resolve the internal structure of the shock. 

Velocity u* corresponds to the maximum value of the curve g{T, u) = in the T-u space, 
as shown in Fig. [71 When the magnitude of this speed is smaller than the magnitude of the 
post-shock speed ui, the temperature increase s monotonically from To to Ti, without presenting 
an isothermal region, for all Prandtl numbers ( Landau &: Lifshitall959l . §88). This occurs in weak 



shocks, called sub-critical, which Mach numbers M2 s smaller than a critical Mach number M2 cr = 
V(37-l)/(3-7)7. 

On the other hand, if the shock is strong, the fixed point diagram looks like the one depicted 
in Fig. [71 In this case, if the Prandtl nui nber is small (as in the s olar corona and our simulation), 



the solution is composed of two regions (jLandau &: Lifshitzjll959l . §88). In the first region, called 



the heat front or thermal conduction front, temperature is increased approximately along the 
g(T,u) = curve, from T2 up to approximately the steady-state value Ti. The second region, 
called isothermal sub-shock or isothermal discontinuity (or jump), consists of an almost isothermal 
{dT/du = 0) change of velocity from uq = Ti/ui to the steady-state value ui. This last region 
corresponds to an isothermal compression of the plasma; most of the compression across the shock 
occurs in this region (see top panel of Fig. [TO]) . The heat front increases the temperature of the 
plasma without changing the density significantly. The compression ratio between the initial density 
P2 and the density at the transition between the sub-shock and the heat front po is bounded between 
the values (7 + l)/2 (limit for a infinitely strong shock) and (87 — l)/(7 + 1) (at the critical Mach 
number), for a Prandtl number equal to zero. For a gas constant equal to 7 = |, this corresponds 
to a density ratio between 1.33 and 1.5. 



5.2. Physical size 

The thickness of the two internal shock regions, sub-shock and heat front, is determined by Eq. 
(j5-5p and (j5-6p . respectively. For the case of strong shocks at low Prandtl number, the thickness 
of the isothermal sub-shock is determined mostly by viscosity, and the heat front length by the 
thermal conduction, as shown below. To calculate these thicknesses analytically, we approximate 
our solution by the case where Pr is zero. In Fig. [H and [5l large and small vertical ticks indicate 
the beginning and end of the thickness of the shock that correspond to the width of the sub-shock 
in addition to the width of the thermal front. 
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In the bottom panel of Fig. [HI the evolution of the plasma-/? of the tube is shown (in the 
reference frame of the shock). The thermal front extends to the right, and the isothermal sub- 
shock extends to the left. The top panel shows the evolution of the length of each region. 



The length (unitless) of the isothermal sub-shock, A/gg^-^_j^Q^ , 
in Eq. ()5-5p . then 
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is calculated replacing T = Ti 



(5-10) 



where integration is calculated from velocity ui to velocity uq = Ti/ui = 2u* — ui (this latter speed 
corresponds to the velocity at the transition between the sub-shock and the heat front as shown in 
Fig. [7]). The plasma-/? at this point can be calculated as /3o = /32U1U2/T2. In the bottom panel of 
Fig. [HI the horizontal dotted line represents /3o. This value is used to approximate the position of 
the transition between the sub-shock and the heat front (zero horizontal coordinate in this Figure). 



The integrand in Eq. (j5-10p diverges at both limits, therefore we will slightly modify them. We 
will define a smaller centered interval corresponding to a 95% of the original integration interval, and 
the integration will be performed with new limits. Defining a small interval Uint = (uo/ui — 1) /40, 
the analytical solution for the length of the sub-shock is prescribed as 
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(5-11) 



where L^°^(T2) = rj{T2) / P2V^\^{T2)Lq is the unit-less mean-free-path of the ions at temperature T2 
(all the lengths are non-dimensionalized to characteristic length Lq); and v^Yii'^^) is the thermal 
speed of the ions at T2. For the simulation presented in Section [4.21 -^^^^'^(72) ~ 0.9 x 10^^. 



Expression (|5-lip presents values for two cases, as shown in the last part of the equation. The 

5 

first one corresponds to temperature independent rj, and the second one corresponds to rj ~ T2. 
The length of the sub-shock is considerably different depending on which model is used for viscosity. 
For the simulation we presented in the last section (temperature dependence of viscosity given by 
Eq. |3-20| ). the ratio between post-shock and pre-shock temperatures was almost one order of 
magnitude, therefore the sub-shock length would be approximately 10^/^ times the one predicted 
by a constant viscosity model. 



For the parameters of our simulation, the bottom part of Eq. (j5-lip gives a length of the 
sub-shock approximately equal to 22 Mm, which is in good agreement with our simulation (see top 
panel of Fig. [H|). In this Figure, the length of the sub-shock (negative positions) increases up to a 
steady-state value well approximated by the analytical solution for Prandtl number equal to zero. 

When a constant viscosity model is used, the length of the sub-shock is smaller than the mean 
free path of the ions for temperature ratios larger than approximately 7 (see Fig. [9]). This is not 
true for the case of more realistic temperature dependencies like the one used in our simulation 
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(r] ~ Ta). For this case, the length of the sub-shock is always more than one order of magnitude 
that of the ion mean-free-path, supporting the use of fluid equations. The length of the sub-shock 
diverges when uq = ui (critical case), as shown in Fig. [9] (mixed dotted and dashed lines), and 
there is no sub-shock for Mach numbers below the critical one. 

To determine the length of the heat front, we will consider the portion of g{T, u) = between 
velocities uq and U2, and replace the corresponding expression for the velocity u = u'* 



in Eq. (j5-6|) . If we define a unitless variable r = — T j (u*)"^, and integral limits T2(M; 
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heat front becomes 



j^^, the integral equation for the (unitless) thickness of the 
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where the "+" sign in the denominator corresponds to super-critical cases, and the "— " sign 
corresponds to the sub-critical cases. In the sub-critical cases the above integral diverges at both 
integral limits, and for super-critical cases, only at T2- The critical case corresponds to ti = 0, 
and the integral is continuous at that point. To define a unique length for all the cases, we will 
integrate expression (j5-12p in a centered interval corresponding to a 95% of the original integration 
interval by defining a small interval Tint = (''"2 — ''"i)/40, analogous to the one use for the isothermal 
sub-shock. 

The case of temperature independent thermal conductivity can be easily integrated since K(r) 
is a constant and can be taken outside the integral. The analytical solution for this case is 
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Here, L^^^'^{T2) = k{T2) / P2V^y^{T2)Lq is the unit-less mean-free-path of the electrons at pre-shock 
temperature T2, and rup and rrie represent the proton and electron masses, respectively. For the 
simulation presented in Section l42l L^^'^^{T2) ~ 0.0021. 

The temperature dependent thermal conductivity case has a slightly more complicated integral 



where k{t) = Kc\u*\ {\ 



For infinitely strong shocks {M2^s = 00)) expression (j5-12p has a 



simpler form since ti(cxd) = and T2(oo) = 1 (the divergence at T2 is removed). In this case the 
integral converges to a finite number, and we see that the length of the heat front is 00 since the 
factor \u*f'/p2U2 (that can be taken outside the integral) increases with the Mach number of the 
shock. The opposite occurs for the constant thermal conduction case. 



For finite shocks, Eq. (j5-12p can be integrated analytically or numerically. For the realistic 
coronal parameters used in our simulation, the length of the heat front is approximately equal to 
93 Mm, in good agreement with our simulation (see top panel of Fig. [8|). This length is of the 
order of the length of a typical coronal loop, and more than one order of magnitude the length of 
the mean-free-path of the electrons at the post-shock temperature. Therefore, for a real coronal 
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situation, the heat front wih probably never reach its steady-state value since it will arrive at the 
end-points before completing its development. At one Alfven time (46 s), the length of the heat 
front is already approximately 28.0 Mm long (the evolution of the heat front is shown in the top 
panel of Figure [8]). 



5.3. Justifying a Fluid Treatment 

In Fig. [21 the lengths of the isothermal sub-shock and heat front for temperature independent 
and temperature dependent models are plotted as function of the post-shock to pre-shock tempera- 
ture ratio. All the lengths are scaled to the electron mean- free-path at temperature T2. Solid lines 
depict the length of the heat front, mixed dotted and dashed lines show the length of the isothermal 
sub-shock, and dashed lines display the particle's mean-free-path at Ti. 

Despite its complicated analytical form, the length of the heat front for the temperature de- 
pendent case follows a simple curve that for relatively strong shocks has a shape similar to the 
electron mean-free-path at Ti (dashed line). This gives some insight into the physics involved. 
The electron mean free path increases with post-shock temperature, as (Coulomb potential), 
while the thermal conduction also increases with post-shock temperature as T2. The interplay 
between particle's mean-free-path and transport coefficient's temperature dependence determines 
the thickness of the shock. Whenever the diffusive coefficient's growth with temperature is faster 
than the particle's mean free-path, the thickness of the shock is going to be larger than the parti- 
cle's mean free path. Heat is transported larger distances ahead of the shock as the temperature 
ratio increases, maintaining a heat front length always larger than the mean free-path. A similar 
analysis can be done for the sub-shock and the ion mean-free-path. When temperature independent 
transport coefficients are used, the results are several orders of magnitude smaller than the corre- 
sponding mean free paths, which in the past created some concern about the use of fluid equations 
for shock studies. 

The validity of assuming fluid equations, as opposed to kinetic theory, can also be a concern 
if the ratio be tween the diffus ion heat flux and the value for streaming particles is larger than a 



certain value (jCampbelllll984l ). We tested the collisionality of the plasma in our simulation, and 
obtained ratios smaller than 0.004. This, and the above paragraph conclusion, support the use of 
continuum fluid equations. 

For the short times involved in a flare, the shocks may not achieve their state state. The 
evolution of the isothermal sub-shock and the thermal front toward this state can be approximated 
by an exponential (see top panel of Fig. [8]). With this approximation, the characteristic time for 
the development of the sub-shock is approximately equal to 532 s, and its initial (maximum) speed 
is less than 10% that of the ion thermal speed. The characteristic time for the heat front is of the 
order of 188 s (four times the Alfven time), and its initial development speed is 3% of the electron 
thermal speed. We have tested the speed of the thermal front for a range of reconnection angles 
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between C = 20 degrees to C = 45 degrees. For all these cases, this speed is between 3% to 5% that 
of the electron thermal speed. 

The electron number density and temperature evolution toward the steady state of our simu- 
lation is shown in Fig. [TOl The electron number density approaches its steady state from above the 
Rankine-Hugoniot value (top panel), and temperature approaches from below (middle panel). A 
short time after reconnection, the two colliding fluids accelerated at the bends get over-compressed 
at the center of the tube (generating densities well above the expected value for infinitely strong 
shocks), and they relax towards the equilibrium value. From the top panel of the Figure, it is 
evident that most of the compression occurs at the sub-shock. The middle panel shows that most 
of the increase in temperature occurs at the thermal front. From the bottom panel of Fig. [HI we 
see that pressure achieves its steady state in a characteristic time much smaller than the ones for 
density and temperature. 

The bottom panel of Fig. [10] shows the volumetric heating rate. Most of the cooling occurs 
at the sub-shock where the plasma is isothermally compressed. Some cooling and heating occur 
at the thermal front. In this region, the plasma increases its temperature without considerably 
changing its density. In the sub-shock region, the velocity changes are large (see Fig. [7|), and it 
would be expected to have a strong positive contribution to the heating, however the cooling and 
heating are dominated by thermal conduction since the Prandtl number is small. The cooling comes 
from the negative curvature of the temperature to the power of 7/2 (not shown in the Figure) as 
described in Eq. ()4-9p . Even for a small temperature gradient in the sub-shock, the heating is 
negative. Without viscosity, however, it is not possible to develop the needed sub-shock to connect 
the pre-shock state to the post-shock state. On the other hand, it is possible to only have viscosity 
and no thermal conduction (limit of Prandtl number equal to infinity). 

6. Conclusions 

We have presented new general fluid equations that describe the dynamics of thin flux tubes. 
These tubes are assumed to be isolated from their surroundings, except by total pressure balance 
with their exterior (their small local radius permits almost instantaneous equilibrium with the 
surroundings via fast-magnetosonic waves). This assumption neglects any drag force due to the 
external fluid. The ideal part of these equations is applicable to a wide range of conditions for 
external magnetic fields and pressure. We also proposed realistic field aligned diffusive terms for 
viscosity and thermal conduction, valid for strongly magnetized plasmas, that include their strong 
dependence on temperature (~ Ta). 

We solve these equations numerically for the reconnection dynamics of a thin flux tube em- 
bedded in an uniform background with small plasma- /3 (/3 = 0.007). This tube was assumed to 
have been reconnected at an impulsive, small patchy region in a current sheet with skewed mag- 
netic fields. Due to magnetic tension, the tube retracts at Alfvenic speeds reducing its length. 
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Realistic coronal parameters (large Reynolds number and small Prandtl number) were used in the 
simulations. 

As the tube evolves, the initially stationary plasma is deflected at the tube bends (rotational 
discontinuities) and accelerated to super-sonic speeds toward the center of the tube where they 
collide. The small, but non-zero thermal pressure at the center of the tube stops the inflows, 
and the two colliding flows generate gas-dynamic shocks that move outward from the center of the 
tube. This shows the importance of including the parallel thermal pressure gradient in the thin flux 
tube equations since there is no magnetic force in that direction. The magnitude of the magnetic 
fleld remains constant across the shocks (these shocks differ from the Petschek case since they are 
not switch-off shocks). The presence of diffusive terms in the thin flux tube equations provides a 
mechanism for the shock development. 

In the simulation, temperature rises almost one order of magnitude from the pre-shock value 
(1 MK) to the post-shock value (~ lOMK) across the gas dynamics shocks. Nevertheless, only a 
small fraction (less than 10 %) of the available magnetic energy is converted to thermal energy. The 
magnetic energy of the tube decreases with time as the tube shortens, and is mostly converted to 
kinetic energy at the bends. This suggest that the rate of the reconnection mechanism has nothing 
to do with the final temperature achieved. This final value only depends on the strength of the 
shock determined by the initial angle between the fleld lines. 

The inner structure of the shock is related to the transport coefficients and strength of the shock 
(Mach number). For solar coronal conditions (strong shocks at low Prandtl numbers), the inner 
structure of the shocks present two sub-regions, a heat front where most of the temperature increase 
occurs (its width determined by thermal conduction), and an isothermal sub-shock where the fiuid 
is compressed (its width depends on the viscosity). For these coronal conditions, the existence of a 
shock transition with thermal conduction alone is not thermodynamically permissible. We estimate 
the thickness of these regions assuming zero Prandtl number and solving the governing differential 
equations analytically. The thicknesses of the sub-shock and heat front are many times the mean- 
free-path of the ions and electrons, respectively. This result, and our calculations regarding the 
free-streaming heat fiux, lead us to conclude that the plasma is sufficiently collisional to use fiuid 
equations. 

The higher density post-shock region should be readily observable. Post-shock density is almost 
four times greater than the initial background density, leading to sixteen times more emission 
measure (emission measure scales as the square of the density). The size of this region increases at 
a rate equal to the speed of the shock (Eq. |5-3j ). and at one Alfven time would be of the order of 



10 Mm. There are several types of observatio i is showing bright features of this scale (iMasuda et al 



1994 



Tsuneta et al.lll997l : IWarren et al.lll999l : I Jiang et al.ll2006l ). The radius of our tube is assumed 



small compared with this length, suggesting that any observable features will consist of multiple 
reconnection events in a small region, and several thin flux tubes being retracted outward from 
that region. 
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The jump condition for thermal pressure is achieved quite rapidly. On the other hand, density 
and temperature achieve their steady-state values on characteristic times longer than that for ther- 
mal pressure. The latter approaches from lower values than the corresponding Rankine-Hugoniot 
post-shock temperature, and the former from higher values. For large Mach numbers (large initial 
reconnection angles), this temporary density overshoot may be larger than the superior limit for 
the Rankine-Hugoniot post-shock density ratio = ^pij)- The value of this ratio may exceed 
one order of magnitude which would manifest as even stronger X-ray emissions near the reconnec- 
tion site, that would decrease in intensity as the tube evolves. Most of the cooling occurs in the 
isothermal sub-shock region, and most of the heating in the thermal front. 

For real coronal situations, the steady-state thickness of the shock might never be achieved 
since the shock will arrive at the end-points before completing its development (the steady-state 
thickness value is of the order of the length of the loop) . For a wide range of initial pre-reconnection 
angles (we tested angles ranging from ( = 20 degrees to C = 45 degrees) the speed at which the 
heat front develops is between 3% to 5% that of the post-shock thermal electron speed. The heat 
front barely compresses the plasma, since most of the compression occurs in the isothermal region. 
The speed at which these fronts move along the legs of the tube depends on the reconnection angle. 
For a magnetic field of Bq = 10 G, and electron number density no = 10^ cm"^, the range in speed 
for the thermal fronts is between 450 km to 650 km s~^ for angles between ^ = 20 degrees to 
C = 45 degrees. For larger angles, it is possible for the heat front to overcome the bends. 

Heat fronts are possibly relevant for chromospheric evaporation in flares. This is not directly 
treatable in our model since it is confined to the current sheet whose edge lays high above the 
chromosphere (see Fig[T]). We expect the heat fronts found in our model to continue along the field 
lines towards their footpoints in the chromosphere. The temperature at the heat front increases 
approximately up to the post-shock Rankine-Hugoniot value, which could be one order of magnitude 
larger than the initial temperature of the loop. At these high temperatures chromospheric material 
can be evaporated, and X-ray emissions from the foot points might be observed. The upward flow 
of evaporated chromospheric material would eventually encounter the sub-shock that is descending 
along the legs of the tube. What would happen after is a subject for further investigation. 

We neglected radiation in our model. The radiative cooling time at the pre-reconnection 
initial conditions is several orders of magnitude larger than the Alfven time. The density increase 
by as much as a factor of four would by itself decrease the radiative cooling time. Temperature, 
however, also increases one order of magnitude which dramatically increases the cooling time since 
the r adiative loss function decreases with i ncreasing tempera ture at temperatures higher than 1 



MK, (ICook et al.lll989l : lMartens et al.ll200d : lRosner et al.lll978l ). Therefore, radiation is a negligible 



effect in the reconnection process. 

Our model assumed a very simple initial geometry: uniform background conditions with pre- 
reconnection field lines that form an angle at both sides of the reconnection . In future work, we will 



investigate a more realistic model with, for example, a Green-Syrovatskii (|Greenlll965l : ISvrovatskii 
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Il97ll ) current sheet with a guide field. This model presents new phenomena as the magnetic field 
magnitude changes along the tube. 

This work was supported by NSF and NASA. We thank Terry Forbes for useful discussions. 

The research in this paper is part of a dissertation by S. Guidoni to be submitted to the 
Graduate School at the Montana State University in partial fulfillment of the requirements for 
completion of a Ph.D. degree. 
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Fig. 1. — Flare current sheet cartoon. The soUd rectangle represents the current sheet plane 
located above the post- flare arcade (gray parabolic loops at the bottom). The black lines depict 
some magnetic field lines on the side of the current sheet that is closer to the viewer, and the dark 
gray ones correspond to field lines on the back side of the current sheet. The lower plane of the box 
represents the solar surface. The small sphere in the current sheet shows a patchy reconnection 
region. Only a small bundle of field lines that intersect this region are going to reconnect. The 
dashed rectangle in the plane of the current sheet corresponds to the region shown in Fig. [3l 
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Fig. 2. — Frenet coordinate system. The axis of a thin magnetic flux tube is represented by a 
central magnetic field line (dashed line) parameterized by vector R(/). The unit vector b is parallel 
to the local magnetic field, unit vector p points in the direction of the curvature vector k, and unit 
vector n completes the system. The two ellipses are the cups of a small tube segment of mass Sm. 
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Fig. 3. — Current sheet cartoon. (a):The two tubes represent the front (sohd contour) and back 
(dashed contour) pre-reconnection tubes. Thick arrows indicate the direction of the magnetic 
field. The sphere in the middle represents the small reconnection region. The tubes are initially 
part of the uniform background equilibrium, (b): Reconnected V-shaped tubes a short time after 
reconnection. The angle between the initial tubes is 2C,. Due to magnetic tension they will retract 
in opposite directions. The end-points of the flux tubes are located near the edges of the current 
sheet. 
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Fig. 4. — Left panel: Simulation results for the temporal evolution of a downward moving recon- 
nected tube. Each horizontal segment corresponds to a different time. Vertical ticks account for 
the thickness of the two outward moving gas-dynamics shocks. The width of the shock evolves up 
to a steady value. Triangles indicate bend positions. Right panel: right bend y-position as function 
of time (same times chosen for the left panel). The solid line correspond to a line with a slope equal 
to the y-projection of the Alfven speed at the reconnection point, vac- Both panels share the same 
y-axis scale 



-34- 




-4000 



4000 



Fig. 5. — Top panel: Dashed line represents initial tube. Solid line represents the tube at time 
equal to 2550 s (vertical scale is compressed). Horizontal arrows delimit regions 1, 2 and 3 (post- 
shock plasma, inflow region, and unperturbed leg of the tube, respectively). Arrows show velocity 
directions in the tube (not to scale). Plasma is deflected in the bisector direction at the bends. 
Large and small vertical ticks indicate the beginning and end of the width of the shock, respectively. 
Bottom panel: the left side depicts temperature profile in half of the tube, and the right side, the 
plasma-beta profile for the same time as in the top panel. Dashed lines show the theoretical steady- 
state Rankine-Hugoniot post-shock values. Crosses are the same as vertical ticks in the top panel. 
In both panels triangles indicate the position of the bends. There is no change in the state variables 
at the bends (Alfven waves). The tube has mirror symmetry in the x-direction. Both panels share 
the same x-axis scale 
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Fig. 6. — Energies' temporal evolution for the simulation presented in Section 14. 2[ Each line 
represents the temporal evolution of the corresponding specific energy integrated over the entire 
tube. The percentage of magnetic energy converted to thermal is small (less than 10 %). The rest 
is converted to kinetic energy. 
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Fig. 7. — Fixed Point Diagram for super-critical Mach number equal to 5.3 (the same one than 
for simulation in Section 14. 2p . Temperature (ordinate) and velocity (abscissa) are scaled to the 
pre-shock values. For this case, \u*\ > \ui\. Top and bottom horizontal dotted lines represent post- 
shock and pre-shock temperatures, respectively. The leftmost vertical dotted line corresponds to the 
post-shock speed, and the rightmost dotted line the pre-shock speed. The dashed line represents the 
curve g{T, u) = 0, and the mixed dotted and dashed line depicts the h{T, u) = curve (the latter is 
the limit solution for an infinite Prandtl number). Each solid line represents a numerical solution 
for equation (|5-9p for different Prandtl numbers (increasing thickness of the line corresponds to 
increasing value of the Prandtl number). The shown Prandtl numbers are 0.01,0.03,0.1,0.3, 1.0, 
and 3.0. The second vertical dotted line from the left shows the characteristic speed u* (maximum 
of the curve g{T,u) = 0), and the second dotted line from the right shows velocity uq (separation 
between sub-shock and heat front for Pr = 0). Asterisks corresponds to grid points inside the shock 
for the simulation presented in Section [4.21 (Pr = 0.01) for time = 2550 s (the same time chosen 
for Fig. [5]). The internal structure of the shock is well resolved by the DEFT program (there are 
several grid points in each shock section). 
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Fig. 8. — Internal shock structure in the reference frame of the shock, for the simulation presented 
in Section 14. 2i The origin is the transition between the sub-shock and the heat front, indicated 
with a solid circle. Both panels share the same x-axis. Bottom panel: plasma-/? as function of 
position in the tube for different times (64 s, 159 s, 319 s, 478 s, 956 s, 1913 s, 2550 s, 2869 s, and 
3188 s). The dotted line represents the earliest time. Horizontal dotted line represents the value 
of Po that corresponds to temperature Ti, velocity uq, and density po- Horizontal mixed dotted 
and dashed line represent the theoretical post-shock values. Plus signs indicate the position of the 
heat front (they correspond to small vertical tick marks in top panel of Fig. [5]). Crosses indicate 
the position of the isothermal sub-shock (they correspond to large tick marks top panel of Fig. [S]). 
Top panel: Temporal evolution of the heat front (solid line) and sub-shock (dashed line) lengths. 
Plus signs and crosses indicate the same positions as in the bottom panel. Leftmost vertical dotted 
line is the theoretical value for the length of the sub-shock calculated from equation (j5-lip . and 
rightmost vertical dotted line is the theoretical value for the length of the heat front. 
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Fig. 9. — Different lengths scaled with the pre-shock electron mean- free-path LF'^^'^iT'i)^ as function 
of the ratio of the post-shock to the pre-shock temperature. The second axis at the bottom shows 
the corresponding Mach number M2,<j. Vertical dotted lines represent cases M2^s = 1 (no shock), 
M2^s = M2^cr (critical case), and M2^s = 5.3 (shock strength for the simulation presented in Section 
4.2p . The lengths in this figure correspond to the same Reynolds and Prandtl number used for the 
simulation presented in Section 14.21 Bottom horizontal dotted line represents the ratio of the ion 
to electron mean free path at T2, and top horizontal dotted line is the unity (electron mean free 
path at T2). Solid lines depict the length of the heat front A/jjp (calculated from Eq. |5-13| and 
|5-12j ). for the cases labeled to the right. Mixed dotted and dashed lines represent the length of 
the isothermal sub-shock A/gg (calculated from Eq. |5-llj ). for the same cases, and dashed lines 
show the mean free path of electrons and ions at post-temperature Ti. 



-39- 




-40 -20 20 40 60 80 100 

(Mm) 



Fig. 10. — Top panel: Electron number density evolution in the shock frame for times t = 63 
s, 159 s, 319 s, 638 s, 1913 s, and 3156 s. Middle panel: Temperature evolution in the shock 
frame for the same times as in the top panel. Bottom panel: Volumetric heating rate in the 
shock frame for the same times as in the top panel. In these three panels, plus signs indicate the 
position of the heat front, crosses indicate the position of the isothermal sub-shock, and dark dots 
indicate the position of the transition between the isothermal sub-shock and the thermal front. The 
dotted lines indicate the earliest time. Horizontal dashed lines are the post-shock Rankine-Hugoniot 
correspondent values. 



